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FLUX-CORRECTED  TRANSPORT  ALGORITHMS 
FOR  TWO-DIMENSIONAL  COMPRESSIBLE  MAGNETOHYDRODYNAMICS 

1.  INTRODUCTION 

Flux-corrected  transport  (FCT)  [1-9]  was  originally  conceived,  and  has  been  developed  over 
the  years,  as  a  method  for  accurately  solving  the  conservation  equations  of  Eulerian  hydrodynamics 
without  violating  the  positivity  of  mass  and  energy,  particularly  near  shocks  and  other  discontinu¬ 
ities.  Thi6  is  achieved  by  adding  to  the  equations  a  strong  numerical  diffusion,  which  guarantees  the 
positivity  of  the  solution,  followed  by  a  compensating  antidiifusion,  which  reduces  the  numerical 
error.  The  crux  of  the  FCT  method  lies  in  limiting  (‘correcting’)  the  antidiffusive  fluxes  before 
they  are  applied,  so  that  no  unphysical  extrema  are  created  in  the  solution.  The  effect  of  this  flux- 
correction  procedure  is  to  provide  as  accurate  a  solution  to  the  original  equation  as  is  consistent 
with  maintaining  positivity  and  monotonicity  everywhere. 

Recently,  a  prescription  was  given  [10]  for  integrating  the  hydromagnetic  equation  of  magne¬ 
tohydrodynamics  (MHD)  using  FCT  techniques.  The  discrete  values  of  the  magnetic-field  compo¬ 
nents  are  placed  at  the  interface  locations  of  the  finite-difference  grid.  In  this  way,  the  divergence- 
free  character  of  the  magnetic  field,  as  expressed  by  the  discrete  integral  form  of  Gauss’s  law,  is 
simply  and  strictly  preserved.  It  was  shown  that  the  monotonicity  constraint  of  FCT,  on  the  other 
hand,  cannot  be  strictly  enforced  without  producing  an  excessively  diffuse  solution  for  the  magnetic 
field.  A  flux-corrector  was  described  which  constructs  the  total  antidiffusive  flux  as  the  sum  of 
corrected  partial  fluxes.  Each  component  of  the  magnetic  field  contributes  a  partial  flux  to  the 
total,  and  its  contribution  is  limited  independently  of  those  due  to  the  other  field  components.  The 
solution  so  obtained  has  remained  well  behaved  even  in  severe  numerical  test  cases. 

In  this  report,  a  pair  of  efficient,  accurate  algorithms  for  integrating  generalized  continuity 
and  hydromagnetic  equations  in  two  spatial  dimensions  will  be  i‘:scr  bed.  The  low-order,  diffusive 
schemes  of  both  the  fluid  and  field  solvers  have  phase  and  amplit<_  errors  that  are  second  order 
in  the  grid  spacing  at  long  wavelengths,  and  are  absolutely  stable  for  Courant  numbers  less  than 
The  errors  in  the  high-order,  dispersive  schemes  sire  fourth  order.  Several  numerical  examples 
will  be  presented  and  discussed. 


ManuscrifX  approved  July  7,  1*>89. 


2.  ALGORITHMS 


Geometry  and  Conservation  Properties 


The  conservation  equations  of  Eulerian  magnetohydrodynamics  take  the  general  form 

^  +  V  •  ( pv )  =  si  +  s2  •  Vs3  +  V  •  s4, 

—  =  Vx(vxB)  +  Vxs5, 


(1) 


where  p  is  a  fluid  variable  (mass,  momentum,  or  energy  density)  being  time-advanced,  v  is  the 
fluid  velocity,  B  is  the  magnetic  field,  and  si,  s2,  53,  84,  and  85  are  source  terms.  In  the  absence  of 
sources,  Eqs.  (1)  reduce  to  the  continuity  and  ideal  hydromagnetic  equations, 


^  +  y.(Pv)  =  0, 

SB  „  . 

—  =  V  x  (v  x  B) . 


(2) 


The  integral  conservation  relations  corresponding  to  Eqs.  (2)  are 

l  =  -  £"  iS' 

^  J  B-dS  =  f  vxBdl, 


(3) 


where  in  the  first  integral  V  is  any  volume  of  fluid  bounded  by  the  closed  surface  5,  and  in  the 
second  5  i6  any  (open)  surface  bounded  by  the  closed  contour  C. 


The  geometry  of  a  finite-difference  representation  of  Eqs.  (1)  is  shown  in  Fig.  1.  It  has  long 
been  recognized  that  the  continuity  equation  can  be  integrated  conservatively  by  evaluating  the 
flux  densities  pv  at  the  interfaces  of  the  spatial  grid  and  using  a  discrete  representation  of  the 
integral  relation  (3a).  For  example,  a  forward  differencing  of  the  time  derivative  and  an  explicit 
treatment  of  the  spatial  derivatives  in  a  two-dimensional  problem  lead  to 

([pc  -  p°]  V)ii  Ar1  =  (p°vxAx)i_^j  -  (p°vxAx)i+ij 
+  (p°»vAy) ij_l  -  (p°vyAy)ij+i , 

where  p°  and  pc  are  the  values  of  the  fluid  variable  before  and  after  the  convection,  respectively,  Ax 
and  Av  are  the  areas  of  interfaces  normal  to  the  x  and  y  directions,  and  At  is  the  time  increment.  In 
a  summation  of  the  results  (4)  over  any  collection  of  cells  in  the  system,  the  contributions  of  fluxes 
evaluated  at  the  common,  internal  interfaces  cancel  pairwise.  The  integral  conservation  relation 
(3a)  in  its  discrete  form  then  holds  for  every  subvolume  of  the  system. 
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Figure  1.  Geometry  of  the  finite-difference  representations  of  the  generalized  continuity  equations  (left 
panel)  and  hydromagnetic  equation  (right  panel)  of  conservative,  Eulerian  magnetohydrodynamics. 


(3a)  in  its  discrete  form  then  holds  for  every  subvolume  of  the  system.  This  result  also  is  true  of 
the  generalized  continuity  equation  (la),  if  only  conservative  sources  s<  are  present. 

A  conservative  integration  of  the  hydromagnetic  equation  is  similarly  effected  [11,12]  by  placing 
the  components  of  the  magnetic  field  at  the  cell  interfaces,  as  shown  in  the  right  panel  of  Fig.  1. 
The  flux  densities  v  x  B  are  evaluated  at  the  cell  edges,  and  a  discrete  representation  of  the  integral 
conservation  relation  (3b)  again  is  used.  If  vx  and  Bx  both  vanish,  discretizing  in  time  as  before 
yields 

m  -  =  ([vvB°t  -  vxB$  Lx).H._h 

-  (K5x  ~  v*Bi\  > 

(5) 

{[Bey  -  By]  Av)ij+i  At-1  =  {[vxBl  -  vyB°x)  Lx)._v+k 

-  {[vxB°v  -  VyB°x ]  Lx).+^j+^  , 

where  B°  and  Bc  are  the  field  values  before  and  after  the  convection,  and  Lx  is  the  length  of  an 
edge  oriented  in  the  z  direction.  Due  to  the  pairwise  cancellation  of  fluxes  at  common  edges,  the 
integral  relation  (3b)  in  its  discrete  form  holds  for  any  surface  in  the  system.  This  result  is  also 
true  of  the  generalized  hydromagnetic  equation  (lb),  whose  conservative  source  term  ss  is  likewise 
evaluated  at  the  cell  edges.  For  the  special  case  of  a  closed  surface,  all  of  the  edge  fluxes  cancel 
and  the  solutions  obey  the  discrete  equivalent  of 

i  B  •  dS  =  0.  (6) 
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Consequently,  if  Gauss’s  law  in  its  discrete  integral  form  is  satisfied  by  the  field  configuration 
initially,  it  i6  satisfied  for  all  time. 

Thus,  the  conservation  properties  of  both  the  generalized  continuity  and  hydromagnetic  equa¬ 
tions  (1)  can  be  simply  and  strictly  imposed  on  their  discrete  solutions  by  a  judicious  choice  of  the 
finite-difference  representation  used  to  solve  them.  A  pair  of  algorithms  exploiting  this  representa¬ 
tion  in  two  spatial  dimensions  have  been  developed,  and  will  be  described  in  the  remainder  of  this 
section. 

Solution  of  Generalized  Continuity  Equations 

The  algorithm  for  continuity  equations  is  the  natural  extension  to  two  dimensions  of  the  one¬ 
dimensional,  low-phase-error  algorithm  ETBFCT  [4,6].  Guirguis  [7]  previously  attempted  such 
a  generalization,  but  his  scheme  suffered  from  excessive  diffusion.  The  key  innovation  introduced 
here  is  the  inclusion  of  nonvanishing  off-diagonal  terms  in  the  antidiffusion  tensor.  Those  terms 
allow  the  phase  and  amplitude  errors  in  the  high-order  solution  to  be  reduced  simultaneously  to 
fourth  order  in  the  mesh  spacing.  If  only  the  diagonal  terms  are  nonvanishing,  at  least  one  of  the 
errors  must  increase  to  second  order,  as  in  Guirguis’  two-dimensional  scheme. 

All  of  these  earlier  algorithms,  and  the  new  one  as  well,  are  one-level  schemes,  whose  errors 
thus  are  first  order  in  the  timestep.  To  achieve  second-order  accuracy  in  time,  a  predictor /corrector 
integration  is  carried  out  using  the  one-level  scheme  for  both  the  half  and  whole  timesteps.  Zalesak 
[5]  instead  used  a  two-level,  leapfrog-trapezoidal  scheme  whose  errors  are  second  order  in  time  and 
fourth  order  in  space,  and  thus  is  comparable  in  accuracy  to  the  new  algorithm  presented  here. 

The  integration  schemes  consist  of  convection,  diffusion,  flux-correction,  and  antidiffusion 
stages.  For  the  continuity  equation,  two  intermediate  solutions  convected  solely  along  x  and  y, 
respectively,  are  needed  in  addition  to  the  solution  convected  in  both  directions,  viz., 

^  =  ^  +  (7) 

Wj  =  WZ  +  *?-»  ~  +  FZ-i  ~  Fh+k' 

where  Vfj  is  the  initial  cell  volume.  The  convective  fluxes  along  i  are  given  by 

=  Pi+ij  A*’ 
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in  which  the  interface  fluid  density  is  calculated  as  the  average 

Pi+iJ  =  2  + 

and  the  velocity  of  the  fluid  with  respect  to  the  moving  grid  is,  in  cartesian  coordinates, 


=  vxi+iJ  -  Af1  -  *?+4)  . 

Here,  x°+ ^  and  x/+^  are  the  initial  and  final  positions,  respectively,  of  the  i  +  |th  interface  along  x. 
The  interface  areas  Act  are  time-  and  space-centered,  and  can  be  chosen  to  preserve  a  uniform 
field  p  under  an  arbitrary  grid  rezoning  with  zero  fluid  flow  [4,6].  The  choice 

=  \  (A Vi  +  Avl)  > 


where 

satisfies  this  requirement.  Analogous  expressions  are  used  to  calculate  the  y  convective  fluxes. 

The  diffusion  stage  includes  the  compression  due  to  grid  motion  and  ensures  the  positivity  of 
the  low-order  solution, 


4vi  -  WZ  +  ff-H  - F?+il  +  Ft,-k  - r?Hi, 

where  V/j  is  the  final  cell  volume.  The  diffusive  fluxes  along  x  are  calculated  from 

=  v***+ii  (p* i  —  P'+ij) 

in  which  the  diffusion  coefficient  is 

.1,1, 

Vxx  i+$j  —  g  +  3fxt+$>» 

the  signed  Courant  number  for  the  fluid  motion  along  x  is 


***+&  ~  2  {^v/j  +  Uxi+ijA*i+^  At' 


and  the  interface  volume  is 


kU,  - 1  iyi + vi,,) . 


(9) 


(10) 


(11) 


5 


Adding  the  contributions  of  the  source  terms  completes  the  calculation  of  the  low-order  solution, 


PyVy^Py^  +  AtV^y 


+  At-  (ACxi+$j  +  Ali-±j)  r  ij  (*>3  i+$j  ~  s3i-lj) 
+  At  j+lj  ■*<  xt+ij  ~  Axi-$j  S4rt-jj 
+  AV  0+  \  S<  V  *J+  AV  ij-  §  S<  V  >j- £  )  » 


(12) 


(13) 


in  which  the  gradient  of  S3  has  been  taken  along  z,  for  definiteness. 

The  antidiffusion  stage,  if  performed  without  flux-correction,  would  yield  the  high-order  so¬ 
lution, 

t>bvu  =  +  fT-a  -  *”+4j  +  fTj-t  -  ni+i- 

The  antidiffusive  fluxes  along  z  are  given  by 

^i+\j  =  Pxxi+\j  (Pi+ij  ~  Pij) 

+  Pxyi+ii  ~  Pi+$}~i)  Vi+\j' 

and  the  antidiffusion  coefficients  are 


(14) 


1 . 

g  g^i+Jj’ 
1 


(15) 


Piyi+^j  —  2€*i+aJ  CV*'+£i" 

The  Courant  number  at  the  z  interface  for  flow  in  the  y  direction  is  conveniently  calculated  as  the 
average 

cv*‘+§i  =  4  (ev*'i-i  (v*+i i-§  €v  «'.?+£  ev»+ij+^)  > 

as  is  the  convected  solution  at  the  cell  edge, 

P*+\i+\  =  4  Pi+ij  P»j+t  +  ^»+lj+l)  • 

The  cross-derivative  terms,  those  with  coefficients  pxy  and  pyx,  make  the  key  contributions  toward 
reducing  the  error  in  the  high-order  scheme. 

In  the  flux-correction  6tage,  the  first  step  is  to  establish  the  allowed  extrema  of  the  final 
solution  from  the  nine-point  comparison 


Pij  —  max  (pi-lj  — l,P«j  — liPi+lj  — l>Pi— lj>  Pij ,  P«+lj,  Pi— 1  j+l,  Pij+\i  Pi+\  j  +  1  )  • 


(16) 
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Here,  pij  may  be  en  to  be  simply  the  low-order  (monotonic)  value 


Pij  —  Pij  * 

corresponding  to  the  ‘strong’  flux-corrector  originally  introduced  by  Boris  and  Book  [1].  Alterna¬ 
tively,  it  may  taken  to  be  any  physically  well-motivated  upper  bound,  e.g., 


pij  =  max  {p°j,p[j) , 

as  suggested  by  Zalesak  [5].  Analogous  expressions  are  used  to  establish  p™in. 

Second,  those  antidiffusive  fluxes  that  are  directed  downstream  with  respect  to  the  local  gra¬ 
dient,  i.e.,  which  act  to  smooth  the  profile  rather  than  steepen  it,  are  cancelled.  Thus,  set 

Fi+H  =  0  ^  *?+*;  {Pi+ij  ~  Pij)  <  0 

and  i?+ij  (p'+2i  -  pli+ij)  <  0  (17) 

or  FT+^{Pij-Pi-ij)<  0. 

Such  fluxes  contribute  to  the  formation  of  dispersive  ripples,  degrading  the  quality  of  the  solution. 
The  third  step  is  the  calculation  of  the  total  antidiffusive  fluxes  into  and  out  of  each  cell, 

=  max  (i^L^o)  -  min  (l?+ii,o)  +  max  (^“-*’°)  -  min  (^“+*,o)  , 

pt~j  =  ma*  (Fi+^°)  ~  “i®  (F“-jj>°)  +  max  (^i+i’0)  "  “i®  (Fij-h'°)  ’ 

and  the  maximum  allowed  fluxes  into  and  out  of  the  cell, 

q*, = wr  -  p.i) v.;. 

05  =  W,  -  list) 

The  ratios  of  these  two  quantities  determine  the  fractions  of  the  incoming  fluxes  that  can  be  applied 
without  causing  an  overshoot  or  undershoot  in  the  final  solution, 

Rtj  ~  min  , 

(20) 

Rij  =  min(l,g-/p-) , 

respectively. 

Fourth,  the  flux-correction  coefficient  for  each  antidiffusive  flux  is  calculated  as  the  minimum 
fraction  which  causes  neither  an  overshoot  in  the  cell  downstream  from  the  flux  nor  an  undershoot 


in  the  cell  upstream, 


C<+ij  ~ 


_  J  min  {Rt+lj,  R~j) ,  if  Fa+ 1 ,•  >  0; 


min  (/?+■,  RT+ij)  i  otherwise. 
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The  fifth  and  last  step  is  to  reduce  the  fluxes  by  these  fractions,  i.e., 


_  ai  ip  a 

Applying  the  corrected  fluxes  to  the  low-order  solution  yields,  finally, 


(22) 


-  API,  +  run  +  K- i  - 


(23) 


the  desired  solution  to  the  generalized  continuity  equation  (la). 

The  positivity  of  the  low-order  solution  pl ,  combined  with  the  restrictions  imposed  by  the  flux- 
corrector,  ensures  the  positivity  of  the  final  solution  p} .  Eq.  (12)  for  pl-  can  be  explicitly  evaluated 
in  terms  of  the  neighboring  values  of  p°  and  cI)V.  It  is  then  easily  established  that  positivity  is 
guaranteed,  in  that  the  coefficient  of  each  of  the  p°’s  is  nonnegative,  for  a  locally  uniform  flow  field 
satisfying  |€XtV|  <  In  the  extreme  case  of  fluid  concentrated  in  and  flowing  outward  from  a  single 
cell,  i.e.,  a  point  blast  wave,  the  limiting  value  which  ensures  positivity  is  considerably  smaller, 
Icx.vl  <  -3)  a  0.15. 

The  accuracy  and  stability  properties  of  the  algorithm  are  analyzed  by  specializing  to  the  case 
of  a  u inform,  cartesian  mesh  and  a  uniform  velocity  field  [3].  The  exact  solution  to  the  continuity 
equation  for  a  uniform  velocity  is  p(x,t)  =  pvx  -  vt,  0).  During  the  time  interval  At,  an  initial 
Fourier  mode 

p(x,0)  =  ,4e’k '*  (24a) 

evolves  to 

p(x,At)  =  p£e,v<x-vAt).  (246) 

The  exact  solution  (24b)  can  be  compared  with  the  numerical  solutions  (12)  and  (13)  obtained 
when  (24a)  is  used  for  p°.  In  the  long-wavelength  limit  k  — ♦  0,  the  phase  and  amplitude  errors  are 
found  to  be  second  order  in  the  mesh  spacing  for  the  low-order  solution  (12),  and  fourth  order  for 
the  high-order  solution  (13).  The  amplification  factor  for  (12)  is  less  than  unity  for  all  wavenumbers 
k,  and  thus  the  algorithm  is  completely  stable,  for  Courant  numbers  |e*,y|  <  %■ 

Solution  of  the  Generalized  Hydromagnetic  Equation 

The  algorithm  for  the  hydromagnetic  equation  is  similar  in  many  respects  to  that  for  continuity 
equations.  Here,  intermediate  low-order  solutions  B '/  and  B‘vv,  to  which  only  the  Bx-  and  By¬ 
dependent  fluxes  contribute,  respectively,  are  needed.  These  solutions  are  not  divergence-free,  but 
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they  are  used  solely  to  correct  their  associated  partial  antidiffusive  fluxes  and  After 

correction,  these  partial  fluxes  are  combined  to  obtain  the  total  corrected  antidiffusive  flux  F“.  As 
discussed  in  greater  detail  elsewhere  [10],  a  flux-corrector  which  limits  the  total  antidiffusive  flux 
F“  with  respect  to  Bx  and  Bv  simultaneously  tends  to  be  very  restrictive,  producing  an  inaccurate, 
diffuse  flnal  solution. 


In  the  convection  stage,  there  results 


,<»is  -  K.+aK.+i,  +  rStj-i  - 
‘Hi  =  AU»i  +  FZii+i  ~ 
=  B°r.+aK  Hi  - 


where  A°  <+^.  and  A°  y + ^  are  the  initial  interface  areas.  The  convective  fluxes  due  to  Bx  and  Bv 
separately  and  jointly  are  given,  respectively,  by 


F'+ii+i  ~  Bxi  uv*+ij+$ 

^•'+4>+i  =  r"*«+ii+i  u*«+£i+£^>  (26) 

rc  _  r»cy  pcx 

«+|j+£  *+ai+^' 

The  edge  magnetic  fields  are  calculated  as  the  averages 

Bxi+$j+±  =  2  (^i  »+£j  +  -®i<+ii+i)  ’ 

5v»+ii+*  =  2  +iBv‘+ii+i)  ’ 

and  the  velocity  of  the  fluid  with  respect  to  the  moving  grid  is,  in  cartesian  coordinates, 

“Ii+3J+J  =  v*»+|i+§  ~  1  (xf+i  ”  1<+i)  ’ 

=  vvi^M~Arl  W+i  “  »J+i)  • 

The  edge  lengths  1  are  to  be  time-  and  space-centered. 

The  diffusion  stage  ensures  the  monotonicity  of  the  low-order  solution  and  includes  the  com¬ 
pression  due  to  grid  motion, 


+  &»- i  -  fSiJH' 

pdx  a/  _  pcx  to  ,  pdy  pdy 

DV'i+\Avij+\  "»«i+iAvu+§  +  «-§j+£  *+£i+§’ 

Biij+iAia+i  -  Bln+iA°vij+i +  F?-a+i  -  F?+a+i- 
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Here,  A^.+  1j  and  A*  ^  are  the  final  interface  areas.  The  diffusive  fluxes  are  given  by 
F*+i>+}  =  ^vvi+iJ+i  ~  Bli+$j+ 1)  A}xi+^.+^, 


Fd\ 


•+ii+i  ^vi+ii+ij  •Av*+ii+i’ 

rd  _  T?dy  rulx 

v+ii+i  «+$j+i  v+£j+£- 


DO  DO 

-  J.  1  ""  -O..  .vi 


\)  Aln 


For  the  diffusion  along  x,  the  coefficient  is 


-  1  1  j 

VXX%+\j+\  -  g  +  gCi+ii+i* 


the  signed  Courant  number  is 


(28) 


(29) 


£*<+*j+i  -  2  l  Af 


»«'i+  i  v  •'+!>+$  . 


u*«+§j+$  ^  <  +ii+i 


and  the  edge  area  is  the  average 

^vi+ii+i  =  2  (^vii+i  +  y4v<+ij+*)  ■ 

Analogous  definitions  apply  along  y.  Adding  the  contributions  of  the  source  term  yields  the  low- 
order  solutions, 

Bxi+\jAii+\i  ~  BtUjAii^  ~  F''+M-l  +  F'+ii+i' 


B»ii+iAlij+}  ~  BiXii+iAlij+$  +  Fi-i)+i  ~  F'+$i+f' 
B*i+h)Azi+$j  =  B*i+$3Aii+$j  ~  F'+±J-$  +  Fi+$j+$' 

KiJ+iAl:Hi  =  BiiJ^AliJ+i  +  F:~^  -  F‘ 


v+£j+£’ 


where 


=  ^*^*t+§j+£5S**+£j+£ 

i6  the  edge  flux  due  to  the  source  term  s5. 

The  high-ord-ir  solutions  are  given  by 

Bh  •  =  B'*»iJAUij  -  Fi+ij-i  +  r?+w 

'''■  hAlij+h  =  Blvij+tAlij+±  ~  F°+ti+f 


(30) 


(31) 
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The  uncorrected  partial  and  total  antidiffusive  fluxes  are 


J-i'tiT 

(^+!j+1  -  *2+i<) 

1  ^xi+ij+i 

+/iyx»+Jj+J 

pay 

•+*;+* 

=Mx*t+Ji+J 

1  RU  dcx  ^ 

ViV+i>+*  ~  BvHH> 

+/*xy  «+!;+§ 

(Bl% ii+.  -  BTi+ii) 

I7»a  _  J7»ay  T-»ax 

•+  ii+i  ’ 


where  the  coefficients  for  antidiffusion  along  x  are 

_  1  1  2 

A*xx  •+$>+$  g  gfxt+$j+£> 

1 

“  2€*,+  aJ+a  €»«+|i+i’ 

The  convected  solution  for  Bx  at  the  y  interfaces  is  calculated  as  the  average 


r>cx  _  f  ( ncx  I  nex  i  tjcx  ,  d cx  'l 

4  Bxi-±j+i  +  ”x»+£j  +  «+ ji+i /  * 


(32) 


(33) 


The  coefficients  for  antidiffusion  along  y  and  the  convected  solution  for  By  at  the  i  interfaces  are 
calculated  analogously. 

In  the  flux-correction  Btage,  the  permitted  extremal  values  of  the  solutions  are  calculated  from 

Bx*+\j  ~  max  (^x«+$j-i»-®x»+$i>-®xt+£j+i)  >  ^ 

Bv*i+\  ~  max  (Bv'-l}+$'Bv*:+i,  Bv>+ij+);')  ■ 

For  the  strong  flux-corrector,  set 

Bxi+$j  — 

=  Byij+^ 

and  for  the  weak  flux-corrector, 

Bxi+hj  =  max  (B°i+ij,B'xxi+i})  , 

Byij+h  =  max 

5"*'."  i  .  and  Bm‘D,  ,  are  calculated  similarly. 

*  *+  5 )  V  'Jt  2 
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Those  partial  antidiffusive  fluxes  that  are  directed  downstream  are  cancelled  according  to 


rpax 

*.-+*;+*  - 

=  0  if 

jpax 

*+£i+i 

(Bxi+$j+l 

and 

jpax 

(B‘xi+$j+2 

or 

rax 

f»+§  i+i 

(^i  - 

JPaV 

.+*;+*  ‘ 

=  0  if 

rav 

*+ij+i 

{Bvi+lj+i 

and 

pay 

•+*;+* 

{?vi+1}+\ 

or 

pay 

KVi  - 

—  R^r 

■D*i+£j+ 1, 


<  o 


(35) 


-  B,v 

V<+lj+* 


i) 


<  0 


Then,  the  ratios  of  the  maximum  partial  antidiffusive  fluxes  into  each  interface, 

K+i)  =  fe-t’1)'”*  >+i'°)  • 

Pt+i 

to  the  maximum  allowed  fluxes  into  the  interfaces, 

<Kni  =  (*%}»- *?»»)<>+»’ 


(36) 


(37) 


(38) 


determine  the  fractions  of  the  incoming  partial  fluxes  that  can  be  applied  without  causing  over¬ 
shoots, 

K+i = (i.e;+j/cj+s)  ■ 

The  fractions  R~+^j  and  outg°ing  fluxes  that  can  be  applied  without  causing  under¬ 

shoots  are  calculated  similarly. 

The  flux-correction  coefficients  for  the  partial  antidiffusive  fluxes  now  are  determined  as  the 
minimum  fractions  which  cause  neither  overshoots  in  the  interfaces  downstream  from  the  fluxes, 
nor  undershoots  in  the  interfaces  upstream, 


CU\i+\ 


Cl, 


j 

fmin  (^+Jj+,,Bi+ij.), 

iff",., 

•+  5J+5 

>  0; 

i 

otherwise; 

_  j 

fmm(a+i>+i,a-+i), 

,+  5r+3 

>0; 

[min  (aJ+J,a-+1J+J), 

,  otherwise. 

(39) 
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Reducing  the  partial  fluxes  by  these  fractions  and  combining  them  to  obtain  the  total  corrected 
antidiffusive  flux  yields 


fpax  _  /~vx  nox 

f'°V  _  /~>V  J?aV 

fo  _  £Say  fa* 

~  *.+*;+*  ~ 

Finally,  the  corrected  fluxes  are  applied  to  the  low-order  solutions  to  get 

~  +  *<+*>+*’ 


(40) 


=  b[.::^aL 


+  F? 


-  F? 


(41) 


'»«+*"»«+*  "  ^vij+hva+i  T  i+i  _  r  <+*;+*’ 

the  desired  solution  to  the  generalized  hydromagnetic  equation  (lb). 

The  accuracy  and  stability  properties  of  this  algorithm  are  identical  to  those  for  the  algorithm 
for  continuity  equations.  A  key  step  in  the  analysis  is  the  elimination  of  B°  in  favor  of  B°,  which 
is  effected  by  using  the  divergence-free  constraint.  For  the  Fourier  mode 


B(x,0)  =  B°keikx, 


the  discrete  form  of  Gauss’s  law  becomes 


sin  0 


*Bo 


6100],  o  _n 
x  k  +  a., 


Ax  x*  '  Ay 

where,  e.g.,  Ax  is  the  mesh  spacing  along  x  and  0X  —  *ArxAx.  There  results 

b1*  B‘.:- 


Es±  =  fsi  = 
b °lr 

V  k 


B°xU 


pV 


so  that  the  errors  in  the  magnetic  field  are  second  order  for  the  low-order  solution  (30)  and  fourth 
order  for  the  high-order  solution  (31),  and  the  algorithm  is  completely  stable  for  Courant  numbers 
less  than 


3.  EXAMPLES 

The  new  FCT  algorithms  described  in  the  previous  section  have  been  programmed  in  single¬ 
precision  FORTRAN  and  optimized  for  use  on  NRL’s  Cray  X-MP/24.  Internal  storage  in  these 
routines  is  minimized  by  stepping  through  the  mesh  one  line  at  a  time,  so  that  no  full,  two- 
dimensional  work  arrays  are  required.  Consequently,  the  modules  use  little  more  memory  than  a 
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MASS  DENSITY 


FCT  FLUID  SOLVER 


Figure  2.  The  mass  density  p  for  a  passive  advection  test  is  displayed  in  its  initial  state  (left)  and  in  the 
final  state  produced  by  the  new  FCT  fluid  solver  (right).  Only  the  50  x  50  subdomain  centered  on  the 
slotted  cylinder  is  shown.  The  axis  of  rotation  is  marked  by  the  solid  dots  in  the  plots  of  the  initial  state. 
The  mass  density  is  contoured  at  10%,  20%,  ...,  90%  of  its  peak  initial  value. 


corresponding  FCT  routine  for  one-dimensional  continuity  equations.  In  hydrodynamical  advection 
tests  on  modest-sized  (100  X  100)  problems,  the  multidimensional  fluid  solver  needed  only  about 
10%  more  time  to  obtain  a  solution  than  did  the  one-dimensional  solver  used  in  operator-split 
mode  (3  fi s  cell-1  timestep-1),  and  its  solution  was  highly  superior  (cf.  [5,6]). 

The  new  algorithms  will  be  applied  here  to  some  kinematical  and  dynamical  problems  of 
hydrodynamics  and  hydromagnetics.  The  first  test  is  the  rigid  rotation  of  a  slotted  cylinder,  as 


calculated  with  the  new  fluid  solver.  Zalesak  [5,6]  introduced  this  test  to  compare  operator-split 
and  fully  multidimensional  solutions  for  an  incompressible  flow  field.  Initial  conditions  for  the 
calculation  are  shown  in  the  left  panel  of  Fig.  2.  The  cylinder  rotates  rigidly  in  the  counterclockwise 
direction  at  an  angular  rate  w,  whence 


vx  =  -uy,  vv  =  +ljx, 


where  the  origin  of  the  cartesian  coordinate  system  is  placed  on  the  axis  of  rotation.  The  calculations 
were  performed  on  a  100  x  100  mesh,  with  the  cylinder  initially  centered  in  the  upper  half  plane 
and  its  radius  set  at  15  zones.  A  slot  of  width  six  zones  runs  through  its  center,  with  a  bridge  of 
the  same  width  left  at  the  top,  joining  the  two  halves  of  the  cylinder.  The  Courant  number  is  0.25, 
so  that  1256  timesteps  are  required  for  one  complete  rotation  in  the  plane. 

The  solution  obtained  by  the  FCT  fluid  solver  is  shown  in  the  right  panel  of  Fig.  2.  The 
quantitative  error  in  the  mass  density,  calculated  from 


( 


VS 


analytic  _  numeric 

Hi) 


analytic 

Hij 


* 


is  25%.  Both  the  Ailing  of  the  slot  and  the  erosion  of  the  bridge  are  limited  to  less  than  10%.  In 
contrast,  an  operator-split  integration  with  the  one-dimensional  routine  Alls  the  slot  and  erodes 
the  bridge  by  more  than  20%.  Qualitatively,  the  multidimensional  solution  compares  very  favorably 
with  Zalesak’s  leapfrog-trapezoidal  results,  as  expected  from  the  similar  error  properties  of  the  two 
algorithms. 

Turning  now  to  the  hydromagnetic  tests,  a  commonly  employed  alternative  to  solving  directly 
for  the  magnetic  Aeld  in  MHD  simulations  is  to  use  a  vector  potential  representation, 


B  =  Vx  A, 


whence  the  generalized  hydromagnetic  equation  (1)  becomes 

6A  _  . 

—  =  vxVxA  +  ss. 


In  two  dimensions,  only  the  component  of  the  potential  along  the  symmetry  direction  is  needed, 
and  this  equation  can  be  reduced  to 


dip 

Hi 


+  V-  Vrp  =  $5, 


(42) 


I 
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where  the  flux  function  ip  is  a  product  of  the  symmetry  component  of  the  potential  and  a  coordinate- 
system-dependent  geometrical  factor.  The  flux  function  is  a  convenient  primitive  variable,  because 
algorithms  for  generalized  continuity  equations  can  be  modified  to  integrate  advection  equations 
such  as  Eq.  (42).  Also,  it  is  a  useful  diagnostic  since  it  cam  be  shown  that 

B  •  Vip  =  0, 

i.e.,  surfaces  of  constant  ip  are  magnetic  surfaces. 

In  the  remaining  examples,  direct  solutions  of  the  hydromagnetic  equation  (1)  are  calculated 
using  the  new  FCT  field  solver.  They  will  be  compared  with  indirect  solutions  obtained  by  inte¬ 
grating  Eq.  (42)  for  the  flux  function  using  a  modified  version  of  the  new  FCT  fluid  solver.  The 
discrete  values  of  ip  are  placed  at  the  cell  edges  (*  +  for  consistency  between  the  two 

calculations.  The  associated  vector  potential  then  is  conservatively  differenced  to  yield  values  of 
the  magnetic-field  components  at  the  cell  interfaces. 

The  second  test  is  the  rigid  rotation  of  a  current-carrying  cylinder,  analogous  to  the  slotted- 
cylinder  test  for  the  fluid  solver.  In  cartesian  coordinates,  the  flux  function  and  the  symmetry 
component  of  the  vector  potential  are  identical,  ip  —  Az.  Initial  conditions  for  the  calculation  are 
shown  in  the  top  panel  of  Fig.  3.  The  potential,  magnetic  field,  and  current  density  within  the 
cylinder  are  given  by 


where  (r,<p,z)  is  a  cylindrical  coordinate  system  centered  on  the  cylinder  of  radius  r0  and  peak  field 
strength  Bo .  A  sheath  of  return  current  at  r  =  ro  neutralizes  the  current  carried  in  the  interior,  so 
that  Ax,  B *,  and  Jx  all  vanish  for  r  >  tq.  The  other  parameters  of  the  calculation  are  the  same  as 
for  the  test  of  the  fluid  solver. 


Figure  3.  [Facing  page]  The  vector  potential  A,  (left)  and  current  density  J,  (right)  for  a  passive  advection 
test  are  displayed  in  their  initial  state  (top)  and  in  the  final  states  produced  by  the  new  FCT  solvers  for 
the  magnetic  field  (middle)  and  vector  potential  (bottom).  Only  the  50  x  50  subdomain  centered  on  the 
cylinder  is  shown.  The  axis  of  rotation  is  marked  by  the  solid  dots  in  the  plots  of  the  initial  state.  The  vector 
potential  is  contoured  at  10%,  20%,  . . .,  90%  of  its  peak  initial  value.  The  current  density  is  contoured  at  ± 
10%,  ±  30%,  . . .,  ±  90%  of  its  initial  interior  value,  with  positive  (negative)  percentages  represented  by  the 
solid  (dashed)  lines. 
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FCT  POTENTIAL  SOLVER 


The  solutions  obtained  by  the  FCT  field  and  potential  solvers  are  shown  in  the  middle  and 
bottom  panels,  respectively,  of  Fig.  3.  The  errors  in  the  vector  potentials  are  very  small  in  both 
cases,  6%  for  the  field  solver  and  4%  for  the  potential  solver.  Although  the  field  solver  clearly 
preserves  the  symmetry  of  the  flux  surfaces  better,  yielding  the  better  qualitative  appearance,  its 
solution  is  somewhat  more  diffuse,  which  accounts  for  its  lower  quantitative  accuracy.  The  errors 
in  the  magnetic  fields  are  essentially  identical  for  the  two  solvers,  at  25%,  and  are  the  same  as 
obtained  in  the  similar  hydrodynamical  test.  For  the  field  solver,  this  error  reflects  the  difficulty 
in  advecting  a  discontinuous  function  (the  magnetic  field  at  r  =  ro)  using  an  Eulerian  difference 
method.  For  the  potential  solver,  on  the  other  hand,  the  error  reflects  the  inaccuracies  in  the 
numerical  derivatives  of  a  smooth  function  (the  vector  potential)  which  has  been  accurately  time- 
advanced  by  a  monotonic  scheme.  The  latter  error  is  further  amplified  in  the  higher  derivatives, 
sufficiently  so  that  the  potential  solver  yields  a  current  density  that  is  less  accurate  than  does  the 
field  solver,  by  90%  to  79%.  While  both  solutions  for  the  current  density  show  local  depressions 
inside  the  cylinder,  the  potential  solver  produces  the  more  severe  fluctuations,  including  a  reversal 
near  the  axis  of  symmetry  whose  amplitude  exceeds  70%.  The  solution  provided  by  the  field  solver 
is  both  qualitatively  and  quantitatively  superior,  although  it  again  suffers  in  the  error  comparison 
because  its  profile  is  not  as  sharp. 

The  third  test  is  the  self-similar  spherical  expansion  of  a  strong  shock  wave  and  trailing 
magnetic  bubble  out  of  the  gravitational  well  of  a  star  [13,14].  In  this  dynamical  problem,  the  self¬ 
similar  character  of  the  expansion  is  maintained  by  a  balance  between  pressure  and  magnetic  forces 
in  the  colatitudinal  direction  and  between  pressure,  magnetic,  gravitational,  and  inertial  forces  in 
the  radial  direction.  The  simplest  expansion  is  an  inertial  flow,  whence 

vr(r,e,t)  =  -t. 

At  time  to,  the  shock  and  its  associated  contact  surface  coincide  at  radius  ro;  at  later  times,  their 
positions  are  given  by 

r.(t)  (t\* 

ro  \to)  ’ 

rc(i)  _(t\ 
ro  Mo/’ 

respectively.  A  magnetic  bubble  is  embedded  behind  the  contact  surface,  and  the  flux  function  rfr 
=  r  sin  BAj,  initially  satisfies 

0(r ,0,  fa )=  /Mr2-r)(r-rj)sin20,  if  r,  <  r  <  r2; 

(0,  otherwise, 
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FLUX  FUNCTION 


FCT  FIELD  SOLVER  FCT  POTENTIAL  SOLVER 


Figure  4.  The  flux  function  rp  =  rain  8A+  for  a  self-similar  spherical  expansion  of  a  magnetic  bubble 
embedded  in  a  stellar  envelope  is  displayed,  in  its  initial  and  exact  fined  states  (top)  and  in  the  final  states 
produced  by  the  FCT  solvers  (bottom).  The  domain  r*o  <  r  <  5ro,  0  <  6  <  ^  is  shown  as  projected  against 
the  sky;  the  shaded  region  is  r  <  r0.  The  flux  function  is  contoured  at  10%,  20%,  . . .,  90%  of  its  peak  initial 
value. 


for  ri  =  0.55  ro  and  r 2  =  0.95  r<).  The  simulations  were  carried  out  for  times  2to  <  t  <  5to,  on  the 
spatial  grid  ro  <  r  <  5ro,  0  <  9  <  x/2.  The  600  timesteps  used  correspond  to  an  average  Courant 
number  of  about  0.25,  and  the  grid  spacing  on  the  100  x  100  mesh  increased  linearly  with  r  and 
was  uniform  in  6.  A  predictor/corrector  integration  method  was  used  to  time-center  the  source 
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terms  and  achieve  second-order  accuracy  in  time.  Finally,  the  gravitational  parameter  GMtl/rl 
was  assigned  the  special  value  2/3,  for  which  both  the  mass  density  and  the  pressure  are  continuous 
across  the  contact  surface,  and  the  flux  constant  ipo  was  chosen  to  yield  a  minimum  plasma  /3  (ratio 
of  plasma  to  magnetic  pressure)  of  unity. 

The  flux  surfaces  at  the  initial  and  final  times  are  shown  in  the  upper  panel  of  Fig.  4.  Shown 
in  the  bottom  panel  are  the  solutions  produced  by  the  FCT  field  and  potential  solvers.  They 
are  quite  similar,  both  qualitatively  and  quantitatively,  although  the  field  solver  yields  notably 
smoother  flux  surfaces.  The  errors  in  the  flux  function  (8%)  and  the  field  components  (20%)  are 
essentially  identical  for  the  two  solvers.  As  in  the  passive  advection  test,  the  potential  solver  yields 
the  less  accurate  current  density,  by  errors  of  113%  to  88%.  Both  solvers  produce  local  current 
reversals  in  this  problem,  however  —  the  potential  solver  at  the  90%  level  and  the  field  solver  at 
the  10%  level.  Contour  plots  of  the  current  densities  are  not  shown  because  they  are  not  very 
instructive.  Numerically  generated  small-scale  structures  dominate  the  current  distribution  within 
the  magnetic  bubble. 

Other  numerical  experiments  have  been  carried  out,  and  merit  brief  mention  here.  First, 
calculations  performed  with  the  bare  low-  and  high-order  schemes  of  the  field  solver  exhibited 
much  larger  errors.  In  the  self-similar  magnetic-bubble  test,  for  example,  each  scheme  produced 
an  error  of  about  60%  in  the  flux  function,  compared  to  8%  when  FCT  is  used  to  interpolate 
between  them.  Second,  a  planar  magnetoacoustic  shock  propagated  over  a  distance  of  60  zones 
broadened  from  an  initial  two-zone  discontinuity  to  a  four-zone  transition  from  pre-  to  post¬ 
shock  conditions.  The  error  in  the  magnetic  field,  evaluated  over  those  four  zones  centered  on 
the  shock  front,  was  less  than  5%  for  the  field  solver.  Third,  the  problem  of  the  expansion  of  a 
cylinder  of  hot  plasma  against  a  uniform  magnetic  field,  for  which  no  analytic  solution  is  known, 
was  simulated  for  a  50:1  pressure  ratio.  In  this  test  for  monotonicity  violations,  the  field  solver 
produced  undershoots  and  overshoots  in  the  solution  of  only  a  few  percent,  despite  the  extreme 
conditions  which  prevailed  early  in  the  calculation.  Finally,  Evans  and  Hawley’s  [12]  monotone 
upwind  scheme  for  the  hydromagnetic  equation  was  applied  to  several  of  these  problems.  Their 
algorithm  is  not  as  accurate  as  the  high-order  scheme  of  the  FCT  solver,  and  consequently  produced 
solutions  with  consistently  larger  errors,  albeit  with  smoother  (more  diffuse)  current  profiles. 
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